function [ dat wnorms wnorm wnorms_nai wnorm_nai weightsout ] = get_voxel_recon_timecourse( S )

% [ dat wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse( S )
%
% MWW 2012

sensor_data=S.sensor_data;
source_recon_results=S.source_recon_results;
chanind=S.chanind;
class_samples_inds=S.class_samples_inds;
voxind=S.voxind; %=first_level_results.mask_indices_in_source_recon(indind)
try, reduce_rank=S.reduce_rank; catch reduce_rank=1; end;

if ~strcmp(source_recon_results.recon_method,'none'), % not sensor space analysis
    NK=length(source_recon_results.BF.inverse.W.MEG);
else
    NK=1;
end;

Ntrials=size(class_samples_inds{1},3);
Ntpts=length(find(source_recon_results.samples2use));

if ~strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
    Nrank=size(source_recon_results.BF.inverse.W.MEG{1}{1},1);
    
    if(reduce_rank>1)
        Nrank=reduce_rank;
    end;
    
else
    Nrank=1;
end;

dat=nan(Ntrials,Ntpts,Nrank);
wnorms=nan(Ntrials,Ntpts,Nrank);
wnorms_nai=nan(Ntrials,Ntpts,Nrank);
wnorm=zeros(NK,Nrank);
wnorm_nai=zeros(NK,Nrank);
weights=zeros(NK,Nrank,length(chanind));

for kk=1:NK,                                
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - just compute for all sensors        
        weights(kk,voxind)=1;
    else
        w=source_recon_results.BF.inverse.W.MEG{kk}{voxind};
        
        if(reduce_rank>1)
            tmp=w*w';   
            [u, ~] = svd(tmp,'econ');                              
            w = u(:,1:reduce_rank)'*w;
        end;
        
        weights(kk,:,:)=w;
    end;
   
    if sum(isnan(squash(weights(kk,:,:))))==0,

        weightskk=permute(weights(kk,:,:),[2 3 1]);
        
        wnorm(kk,:)=trace(weightskk*source_recon_results.BF.features.C.MEG{kk}*weightskk');
        wnorm_nai(kk,:)=trace(weightskk*weightskk');

    else,
        disp('NAN weights');
    end;
    
end;


for kk=1:NK,

    weightskk=permute(weights(kk,:,:),[2 3 1]);            
    weightsout{kk}=weightskk;

    for tri=1:Ntrials, % indexes trials

        timeinds=find(class_samples_inds{kk}(1,:, tri)); % time indices for class kk
        if ~isempty(timeinds),
            dat(tri,timeinds,:)=(weightskk*sensor_data(:,timeinds,tri))';
            %wnorms(tri,timeinds)=wnorm(kk);    
            wnorms(tri,timeinds,:)=permute(repmat(wnorm(kk,:),[1,1,length(timeinds)]),[1 3 2]);
            wnorms_nai(tri,timeinds,:)=permute(repmat(wnorm_nai(kk,:),[1,1,length(timeinds)]),[1 3 2]);
        end;
    end;
end;


end

